Analytic Approximations for the Velocity of Field- Driven Ising 

Interfaces 

Per Arne Rikvold 1 and M. Kolesik 2 ' 3 
February 1, 2008 

Abstract 

We present analytic approximations for the field, temperature, and orientation dependences of the 
interface velocity in a two-dimensional kinetic Ising model in a nonzero field. The model, which has 
nonconserved order parameter, is useful for ferromagnets, ferroelectrics, and other systems undergoing 
order-disorder phase transformations driven by a bulk free-energy difference. The Solid-on-Solid (SOS) 
approximation for the microscopic surface structure is used to estimate mean spin-class populations, from 
which the mean interface velocity can be obtained for any specific single-spin-flip dynamic. This linear- 
response approximation remains accurate for higher temperatures than the single-step and polynuclear 
growth models, while it reduces to these in the appropriate low-temperature limits. The equilibrium SOS 
approximation is generalized by mean-field arguments to obtain field-dependent spin-class populations for 
moving interfaces, and thereby a nonlinear-response approximation for the velocity. The analytic results 
for the interface velocity and the spin-class populations are compared with Monte Carlo simulations. 
Excellent agreement is found in a wide range of field, temperature, and interface orientation. 

KEY WORDS: Kinetic Ising model, Solid-on-Solid (SOS) approximation, Microscopic interface structure, 
Surface anisotropy, Surface growth, Interface dynamics, Linear response, Nonlinear response, Monte Carlo 
simulation. 



J. Stat. Phys. 100, in press. |cond-mat/99Q918& 



1 Center for Materials Research and Technology, School of Computational Science and Information Technology, and 
Department of Physics, Florida State University, Tallahassee, FL 32306-4350. 

E-mail: rikvold@csit.fsu.edu 

2 Institute of Physics, Slovak Academy of Sciences, Bratislava, Slovak Republic 

3 Department of Mathematics, University of Arizona, Tucson, AZ 85721. 
E-mail: kolesik@acms.arizona.edu 



1 



1 Introduction 



The appearance of the world around us through sight and touch is largely determined by interfaces between 
phases with different optical and mechanical properties, and the materials properties of multiphase media 
are strongly influenced by the geometry of interfaces that separate their constituent materials jjj], §]. Since, 
in general, the morphology of an interface is determined by the growth process by which it is formed, the 
large-scale structures of growing interfaces have inspired an enormous amount of work in the last few decades 



In comparison to the vigorous interest in large-scale structure, much less attention has been paid to 
intcrfacial structure on a microscopic scale. This is somewhat surprising since the microscopic structure 
limits the interfacial propagation velocity under an external driving force, such as the applied field for a 
magnetic or dielectric domain wall or the supersaturation or supercooling for a crystal surface. It is also 
important for properties such as chemical reactivity and catalytic activity. 

In this paper we consider how the microscopic interface structure determines the growth velocity of 
a simple model surface in a system with nonconserved order parameter: the interface between domains 
of positive and negative magnetization in a square-lattice kinetic Ising ferromagnet with nearest-neighbor 
interactions, which is driven by a field favoring one of the two spin orientations J5|, ^, This model is 
applicable to the kinetics of phase transformation in many magnetic and ferroelectric systems and other 
order-disorder transitions whose kinetics are not inhibited by coupling to a conserved field. It belongs to the 
dynamic universality class of the Kardar-Parisi-Zhang (KPZ) model fa, and the velocity of a macroscopically 
plane interface is expected to be linear in an asymptotically weak field, as is also predicted by the Lifshitz- 
Allcn-Cahn theory || [n|. However, neither theory gives the explicit field dependence, which should contain 
both the average interface orientation and the specific dynamic. Here we derive analytic, approximate 
expressions for the mean velocity as a function of field, temperature, and interface orientation. Our approach 
is based on the concept of spin classes used in rejection-free Monte Carlo (MC) algorithms [l^, p"3| ], 
together with the Burton-Cabrera-Frank Solid-on-Solid (SOS) approximation for the structure of a stationary 
interface |l4|, ^5). While the theory should become exact for asymptotically small temperatures and fields, 
our main purpose is to explore its applicability outside this limited regime. 

The remainder of this paper is organized as follows. In Sec. § we introduce the kinetic Ising model and 
the concept of spin classes. In Sec. [| we summarize relevant aspects of the SOS approximation for the 
structure of a flat, equilibrium Ising interface between two bulk phases of opposite magnetization, which 
we use to obtain analytic approximations for the mean spin-class populations in zero field. These provide 
a linear-response approximation for the velocity of a driven interface. In Sec. ^ we develop an extension of 
the SOS approximation to obtain field-dependent spin-class populations for flat, moving interfaces as well. 
This leads to a nonlinear-response approximation for the velocity. In Sec. ^| the theoretical results for the 
interface velocity and spin-class populations from Sees. || and |J are compared with MC simulations. The 
nonlinear-response approximation gives remarkable agreement with the simulations in a wide range of field, 
temperature, and interface orientation. Section |6| contains a discussion, conclusions, and some suggestions 
for future work. 

2 Model and Dynamics 

The anisotropic square-lattice Ising ferromagnet with nearest-neighbor interactions is defined by the 
Hamiltonian 



Here s x>y — ±1, y runs over a U lattice sites, and H is the applied field. The lattice constant is taken as 
our unit of length. An interface is introduced by fixing s x y = +1 and —1 for large negative and positive 
y, respectively. For concreteness we assume that H > 0, such that the interface moves in the positive y 




(1) 
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direction under an applied field. The implementation of these boundary conditions in our MC simulations 



is discussed in Sec. 5.f 



Approach to equilibrium is ensured by a single-spin-flip (nonconservative) dynamic which satisfies detailed 
balance, such as the Metropolis or Glauber algorithms pl| . Any such algorithm is defined by a transition 
probability, W[s x>y — > — s x , y ] — W[/3AE], where f3 is the inverse of the temperature T (we use units in which 
Boltzmann's constant is unity), and AE is the energy change that would occur if the proposed spin flip 
were accepted. Since there are only a finite number of different values of AE, the spins can be divided into 
classes || |Tf[ , labeled by the spin value s and the number of broken bonds between the spin and its nearest 
neighbors in the x- and y-direction, j and fc, respectively. The spin classes, denoted jks with j, k £ {0, 1, 2}, 
are listed in Table 1 together with the corresponding energies, E(jks), and energy changes, AE(jks). For 
the anisotropic model defined by Eq. (|l|) there are f8 classes In the isotropic case, J x = J y , this reduces 
to 10 classes distinguished by s and the total number of broken bonds, j+k [p"l|| . 

For concreteness and comparison with numerical simulations we here choose the discrete-time Glauber 
dynamic, defined by the transition probability 

W G [s XtV -» -s x . y ] = 1 + e ^ fJAE • (2) 

The Glauber dynamic is mathematically convenient in that the transition probability is a continuously 
differentiable function of AE. In its continuum-time version it has also been shown to correspond to a 
quantum-mechanical S=l/2 system weakly coupled to a large thermal fermion bath Jl7| . However, the 
spin-class populations can be used to estimate propagation velocities with any single-spin-flip dynamic that 
satisfies detailed balance. Time is measured in units of MC steps per spin (MCSS). 

To prevent nucleation of droplets of the stable phase in front of the moving interface [|| , we modify the 
Glauber dynamic by setting the transition rate for any spin which is parallel to all its neighbors (i.e., class 
00s) equal to zero 0, [L8[ fj|). This suppresses thermal fluctuations in both the bulk phases, while the local 



interface structure is reasonably preserved. For moderate fields the interface velocity with this modified 
Glauber dynamic is only slightly less than that obtained with the full Glauber dynamic Jl8|, [H^. In the 
strong-field regime, where the size of a critical droplet of the equilibrium phase is reduced to on the order 
of the lattice constant, the kinetic Ising interface loses its integrity. A conservative analytic approximation 
for this crossover field (often called "the mean- field spinodal" [^0[ |l|) is (for isotropic interactions) |22| 
Hmfsp(T) w a(T)/m cq (T), where cr(T) and m e q(T) are the equilibrium surface tension in the ^-direction 
and the equilibrium magnetization, respectively. In the strong-field regime the SOS approximation and the 
dynamic used here should be considered as a nonequilibrium cluster growth model in its own right. For 
\H\ < Hmfsp{T) they constitute a good approximation for the kinetic Ising model with the full Glauber 
dynamic pi 111. 



3 The SOS Approximation 

The separation of spins into classes forms the basis of several rejection-free MC algorithms [0, [l^, [13], |2^, 
p4| , p5[ . In such algorithms the spin-class populations, n(jks), are continually monitored throughout the 
simulation. Given this information and the transition probabilities of the particular dynamic used, one can 
then calculate the time increments between MC updates. Here we instead obtain analytic approximations 
for the mean spin-class populations for a driven interface moving at a constant velocity, based on the Burton- 
Cabrera- Frank SOS model of the equilibrium interface [O. These populations are then used together with 
the transition probabilities to obtain the mean interface velocity. 

The SOS approximation describes the interface as a single- valued function y(x). For the square lattice 
considered here, the interface is a series of integer- valued steps of height 5(x) parallel to the j/-axis, as shown 
in Fig. [|. The heights of the individual steps are assumed to be statistically independent and identically 
distributed. The probability density function (pdf) is given by the interaction energy corresponding to the 
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\S(x)\ broken J^-bonds between spins in the columns centered at (x — 1/2) and (x + 1/2) as 

p B [5(x)] = Z- X X\ 5 ^ e^ 6 ^ , (3) 

where we have introduced the shorthand X = e~ 2 P J: *. Here 7(0) is a Lagrange multiplier which maintains 
the mean step height at an x-independent value, (S(x)} = tan0, where <fi is the overall angle between the 
interface and the x-axis. The partition function is 

+ °° l y2 

6 ~ 1-2X cosh 7 (0) + X 2 ■ lj 

S— — 00 

Using Z x as a moment generating function for 5 (x), it is straightforward to obtain the explicit expression 

w = (l + X 2 )tan0 + fl 

2X(l + tan0) ' V 7 



where R 
of 0: 



T-X 2 ) tan 2 + 4J*P 



1/2 

. Combining Eqs. and (g) one obtains Z x explicitly as a function 

(1 - X 2 ) (1 - tan 2 0) 

Zx{4>) = ~ — ■ (6) 

x\v>) l + X 2 -R W 

For = this simplifies to Z x (0) = (1 + X)/(l — X). The SOS approximation ignores overhangs and 
bubbles. It is therefore rather remarkable that the surface tension in this approximation, calculated as 
csos = I cos (f>\ [2 J y — Tin Z x {&) + T^f(4>) tan0], yields the exact result for (j) — J2(| and an excellent 
approximation for \<f>\ < ir/4 |27|| . For larger \(j>\ it is more reasonable to use an SOS approximation with 
steps parallel to the x-axis. 

While Eq. (|^) is equivalent to Eq. (72) of Ref. |fL4j| , the implicit form given by Eqs. (||) and (||) is more 
convenient for our purpose of obtaining mean spin-class populations, (n(jks)}. The spin classes compatible 
with this approximation are illustrated in Fig. [p. The mean populations are all obtained from the joint 
pdf for 5(x) and 5(x+l). Since the individual step heights are statistically independent, this is the product 
p x [<S(:r)] -Px+i [^( :c + 1 )]- The symmetry of p x [(5(a;)] under the transformation (x, <fi, 5) — > (—x, —(f), —S) ensures 
that (n(jk— )) = (n(jk+)) for all j and k. (On the right-hand sides of Eqs. (0) - ( |1C| ) below, we have chosen 
s = —1 with the interface oriented as shown in Fig. [j] for concreteness.) 

The SOS picture implies that there is exactly one broken J y bond per unit length in the x-direction, so 
that (n(01s)) + (n(lls)) + (n(21s)) = 1. The calculations of the individual populations are straightforward 
but somewhat tedious, especially for nonzero cj>. In Eqs. (0)-(|^) below we therefore just give the starting 
point of the calculation for each class in terms of p x [5(x)] and the cumulative probabilities P[5(x) < n] = 
^^ = _ 00 p x (<5) and P[S(x) > n] = 1 — P[S(x) < (n — 1)]. The final results are listed in Table 2, both for 
general cf> and for = 0. 

(n(01s)) = P[S(x) > 0] ■ P[S(x + 1) < 0] (7) 
(n(lls)) = P[8(x) < —1] • P[8(x + 1) < 0] + P[S(x) > 0] ■ P[S(x + 1) > 1] (8) 
(n(21s)) = P[S(x) < -1} ■ P[5(x + 1) > 1] (9) 

Flipping a spin in either of these classes (marked * in Table 2) preserves the SOS configuration. 

To obtain the mean populations for classes of spins that are connected to the interface only through one 
(10s) or two (20s) broken J x bonds is more tedious. We found it most convenient first to calculate the joint 
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pdf for n = ra(10s) + n(20s) and n(20s): 

' P[8(x) >-!]■ P[8(x + 1) < 1] for n = n(20s) = 



p[n, n(20s)] = < 



p x [-(n + l)]-P[5(x+l) < 1] 
+P{S(x) > -1] - Px+1 [n+l] 



for n > 1, n(20s) = 



^[-(n + l)] •p a+1 [n(20s) + l] 

+p s [-(n(20s) + 1)] ■p x+ x[n+ 1] for n > 1, 1 < n(20s) < n 



(10) 



P^[-(«+ 1)] ■p a; +i[n + 1] 



for n(20s) = n > 1 



The resulting expressions for (n(10s)) and (n(20s)) are marked f in Table 2. Flipping a spin in one of these 
classes results in the breaking of two Jy-bonds, and thus in the creation of an overhang or a bubble. The 
classes, 02s, 12s, and 22s, are not populated in the SOS approximation, while 00s represents bulk sites which 
have zero transition rates with this dynamic. 

Each column of the interface advances by one lattice constant in the y-direction whenever a spin flips from 
— 1 to +1, regardless of its y-coordinate. Conversely, the interface recedes by one lattice constant whenever 
a spin flips from +1 to —1. The energy changes corresponding to a flip are given in the third column in 
Table 1. Since the spin-class populations on both sides of the interface are equal in this approximation, the 
contribution from sites in the classes jk— and jk+ to the mean velocity in the y-direction is the difference 
between the transition probabilities for spin flips leading to advance and recession: 



(v y (jk)) = W [pAE(jk-)} - W [f3AE(jk+)} 



(11) 



The results corresponding to the Glauber transition probabilities from Eq. (|2|) are given in the last column 
in Table 2. (It is of course trivial to generalize to (n(jk—)) ^ (n(jk+)), but this will not be needed here.) 
The mean propagation velocity perpendicular to the interface becomes 



(v±(T,H,<l>)) 



4,\J2(n(jks))(v y (jk)) . 

3,k 



(12) 



Including the SOS-violating moves, jk £ {10,20}, in principle allows the propagation velocity to exceed 
unity as it becomes possible to flip several spins along a high vertical step. However, such large velocities 
are only observed for the strongest fields investigated here. Restricting the sum to only jk £ {01, 11,21}, 
on the other hand, yields an approximation which excludes SOS-violating transitions and would limit the 
velocity to below unity. 

While the general result for the velocity is rather cumbersome if written out in detail, the special case of 
<fi = leads to a relatively compact formula: 



(vx(T,H,0)) = 



tanh(/3if) 



x 2 



l-X 2 



2X 



1 + X 2 



sinh(2/3J ;c 



cosh(,3H) 



X 2 



1 



sinh(2j3(J y -J :c )) 
cosh(l3H) 



1 



2(1 + 2X) 



sinh(2/3J H ) 
cosh(/3ff) 



(13) 



Here the first line corresponds to transitions that preserve the SOS structure of the interface, while the 
second line corresponds to transitions that create overhangs or bubbles. Comparison with simulation data 
indicate that excluding the SOS-violating transitions leads to significant underestimation of the propagation 
velocity, even for quite moderate fields. This effect is shown in Fig. Band discussed in more detail in Sec. 5.2. 
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Equation (|T^) with J x — J y was presented without detailed derivation in Ref. |T^|. In that work 
the average interface velocity in a kinetic Ising model undergoing a field-driven phase transformation was 
estimated directly from the time evolution of the magnetization after a sudden field reversal. This estimate 
was found to be consistent with Eq. (|i"3|). 

As T — > 0, X — > 0. In this limit, the dominant term in Eq. ( p^ ) is the one proportional to (n(lls)), 
which is simply the density of kinks on the surface. Combining the appropriate entries in Table 2 with 
Eq. (g) and ignoring X everywhere except where X 2 occurs in an additive combination with tan 2 cf>, we get 
the angle-dependent interface velocity for very low T: 



(v x (T^0,H,<f>)) = cos <\> v u " ; T " 2 M v tanh (f3H) (14) 



Vtan 2 4> + AX 2 - tan 2 , 
1 - tan 2 

1 tanh (J3H) for In no » A 



2 | sin 0| + | cos <f>\ 



(15) 



v/tan 2 + AX 2 tanh (/?#) for tan </>< X 



The first line of Eq. (|l5|) corresponds to the single-step model ^ |[ ^9|. In Ref. J(| the same result 
was obtained from a Green-Kubo-like linear-response formula for the interface mobility. In the second line 
we retain only those terms in tan 2 <fi, which dominate in the limit |tan</>| «!,!-> 0. It corresponds to 
the polynuclear growth (PNG) model || ||(| |TJ. Equation (14) provides a correct interpolation between 



the PNG and single-step results as |tan</>| increases from |tan^| <C X to | tan <^>| ^> X. Only the factor 
tanh(/3_Jf), which corresponds to the average velocity of a single step along the surface, depends explicitly 
on the specific dynamic. At higher temperatures than those for which Eq. (|14j) holds, the other spin classes 
contribute to the interface velocity as well, and our approximation goes beyond the single-step and PNG 
approximations . 



4 Nonlinear-Response SOS Approximation 

The velocity estimates obtained in Sec. g were derived from the equilibrium interface fluctuations at H = 
and thus constitute a linear-response approximation. This is satisfactory for sufficiently weak fields. For 
stronger fields, however, the structure of the moving interface is, in general, modified, leading to additional 
field dependences in the velocity. As will be shown in Sec. [| this effect can be significant. In this section we 
therefore develop a mean-field theory for the step-height pdf for a moving, flat SOS interface in a nonzero 
field. 

The step-height pdf for a stationary SOS interface in a nonzero field is E| 

p x [6(x)} = Z-'X^ e bW-^Hx]8 {x ) (16) 

The term containing H adds an cc-dependence to the Lagrange multiplier which determines (<5(x)).Thc 
corresponding ^-dependence in the partition function is obtained by replacing 7(36) by ["f(4>) — 2[3Hx\ in 
Eq. (0), where now tan0 = {5(G)). The geometric structure described by Eq. ( p"6|) is a macroscopically 
curved interface which bulges in the direction of the metastable phase region. For the case of conserved 
order parameter or nonconserved order parameter with an interface pinned at two points, the stationary 
configuration is an equilibrium one. For nonconserved order parameter without pinning, it corresponds to a 
critical droplet of the stable phase. If the average curvature is changed from the x-dependent form given by 
Eq. (|l6[), the interface will move. 

It is well known that the macroscopic, stationary distribution for flat, moving interfaces in the 
KPZ universality class is Gaussian corresponding to a random walk with independent increments. 

Nevertheless, the step heights in several discrete models in this class are known to be correlated at short 
distances |52|, [53| . Here we develop a mean- field theory for the single-step-height pdf of the moving interface 
by ignoring these short-range correlations. Thus we assume that the step heights are statistically independent 
and identically distributed, just as they are for H = 0. In this approximation, the single-step pdf of a moving 
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interface parallel to the x-axis is given by Eqs. (||) and ([!]) with 7 = and an ^/-dependent generalization 
of the parameter X. We now construct a self-consistency equation to determine this parameter, X(T, H). 

For simplicity we consider the step at x = 0, and we take 8(0) > 0. The equations obtained by considering 
other values of x and 8(x) < are identical to the one we derive below. The total transition probability for 
the height of a single step to change from 8 to 8 ± 1 is W[8 — ► 8 ± 1], Relating the single-step transition 
probability to the single-step-height pdf by detailed balance, we have: 

x(t m - Mm + 1] - w[m - m + 1] (m 

K ' ] ~ Po [S(0)] W[S(0) + 1^8(0)] • 1 ' 

To find these transition probabilities, we refer to Fig. ||. In order for 8(0) to increase by one, either the 
spin at x = +1/2 can flip from —1 to +1, decreasing <S(+1) by one, or the spin at x = —1/2 can flip from 
+ 1 to —1, decreasing S(—l) by one. Each of these possibilities is attempted with probability 1/2. For each 
spin-flip direction, the resulting energy change can have two different values, depending on the height of the 
corresponding neighboring step. Analogous arguments hold for the reverse transitions, 6(0) + 1 — > 6(0). The 
energy changes and the conditions on the neighboring step heights are given next to the arrows denoting the 
directions of the transitions in Fig. ^. Expressing the single-spin transition rate corresponding to an energy 
change AE as PF[Ai5] and reading the energy changes and conditions off from Fig. ||, we get: 

W[8(0) -f 8(0) + 1] = -{w[-2H]P[8(+l) > +1] + W[-2H + 4J x ]P[5(+l) < 0] 

+W[+2H]P[8(-1) > +1] + W[+2H + 4J X ]P[8(-1) < 0]} 

W[S(0) + 1 -»■ 8(0)] = ^{W[+2H]P[8(+1)>0]+W[+2H-4J X ]P{S(+1)<-1} 

+W[-2H]P[8(-l) > 0] + W[-2H - AJ x ]P[8(-l) < -1]} 

(18) 

Consistent with the mean-field approximation we calculate the cumulative probabilities with the same 
stationary single-step pdf, obtaining P[8 > 0] = P[8 < 0] = 1/[1 + X(T,H)} and P[8 > +1] = P[8 < 
-1] = X(T, H)/[l + X(T, H)}. The resulting self-consistency equation for X(T, H) is 

= X(T, H) {W[-2H] + W[+2H}} + W[-2H + 4J X ] + W[+2H - 4J X ] 
[ ' ' W[-2H] + W[+2H] + X(T, H) {W[-2H - 4 J x ] + W[+2H - 4J X }} ' 1 1 

This is solved to yield 

X(T,H) = e"^- K 2 ^- 2g -^^ , (20) 

y J 1 W[-2H - 4J X ] + W[+2H -4J X ] J y ' 

where we have also used the detailed-balance relation for the single-spin transition probabilities, VF[+A£']/FF[— AE] 
e- f3 ^ E , to eliminate W[+2H + 4J X ] and W[-2H + 4J X ]. 

The approach described above is also applicable to curved interfaces, and the resulting self-consistency 
equations are cubic in X(T,H). However, except for the stationary curved interface described by Eq. jilf), 
the shapes of such interfaces are not stationary and much more difficult to investigate by simulations. We 
hope to return to these more complicated problems in the future. 

Equation (^0|) shows that X(T, H) depends on the specific dynamic, except for H = 0, where it reduces 
to its equilibrium value, X(T, 0) = e~ 2l3J:c , for all dynamics that satisfy detailed balance. 

Using the Glauber dynamic defined by Eq. (0), we get explicitly 



2(3J - COSh(2/?#) 4- e-WJ* 



e" H " m costirzpu + e 



x g( t, H) = \ : ,,r,:T. ) ;: //l ; ( ^, \ ■ (^d 



e -20j x C osh(2/3i/) 



1/2 
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The Metropolis dynamic, Wm[A£] = min [l, e _,3AB ] , yields 



-2/3j x {cosh(2/3 J ff)} 1/2 for H < 2J X 

\T/.// { i/2 ■»■> 



{ ^ ( "tr;i } foriI>2J, 



Numerically, Xg and Xm are not very different, and they both approach unity from below as H increases 
beyond 2J X . They are shown together vs H in Fig. || 

The approximation used in Ref. @, X{T,H) = e" 2/3J - cosh(/3iJ), can be obtained by "brutally 
decoupling" the steps through fictitiously splitting each spin in two and flipping only half of a spin to 
change the height of a single step. For such a process Eq. ( |l7| ) directly yields X(T,H) — {W[— H + 2J X ] + 
W[+H + 2J X ]}/{W[+H - 2J X ] + W[-H - 2JJ}, from which the result in Ref. @] is obtained for H < 2J X 
using the Metropolis dynamic. It is surprisingly close to the proper mean- field results and is included in 
Fig. U for comparison. 

In their study of a lattice-gas model for three-dimensional crystal growth |34| , Kotrla and Levi introduced 
a single-site dynamic which in Ising language can be described (up to an overall rate factor e 2 " H ) as the 
product of Metropolis for field effects and Glauber for interaction effects: Wkl[—2H — 4JJ = 4JJ 
and Wkl[+2H — AJ X ] = er 2 ^ H Wo,\—\J x \. Inserting these transition probabilities in Eq. (fj()|), one finds 
that they lead to an ff-independent X. This is indeed the case for any dynamic in which the field and 
interaction effects are statistically independent and obey detailed balance separately. It shows that the local 
interface structure of driven interfaces can vary strongly as the dynamics are changed, even within the class 
of nonconservative single-spin-flip dynamics. Dynamics that factorize in this way are known as "soft," in 
contrast to the "hard" dynamics which do not factorize, such as Metropolis and the Glauber rate used here. 
Soft and hard dynamics have been shown to lead to differences in the steady states in a number of other 
nonequilibrium systems as well p5| . 

All the results for the spin-class populations of the zero-field equilibrium interface, which were derived 
in Sec. |^ and are listed in Table 2, can now be applied to obtain a nonlinear-response approximation for 
the propagation velocity of flat, driven interfaces. This simply requires replacing the zero-field X = e ~ 2 P J * 
used in the linear-response approximation by the field-dependent X(T, H) 7 obtained from Eq. (|2(]) with the 
transition probabilities corresponding to the particular dynamic used. For most hard dynamics, the net 
effect is to increase the mean step height, (|<5|), and thus the interface velocity. 

In the next section we show that this nonlinear-response approximation gives good agreement with MC 
simulations of driven, flat Ising interfaces for a wide range of fields and temperatures. 



5 Comparison with Monte Carlo Simulations 

We have compared the analytical estimates of propagation velocities and spin-class populations developed 
above with MC simulations of the same model for J x = J y = J. 

5.1 Simulation details 

To minimize the finite-size effects (see below) , large simulation lattices are needed to accommodate a sufficient 
length of the interface and provide enough room for it to move unimpeded. To achieve long simulation runs it 
is necessary to employ a co-moving simulation box to prevent the interface from hitting the system boundary 
in the y-direction. We therefore used an "active-zone" algorithm which relies on the fact that spins with 
no broken bonds have zero transition probability. The algorithm uses a lattice of size L x in the x-direction, 
and it keeps a list of all spins which have at least one broken bond and thus can flip in the next simulation 
step. Once a spin loses its last broken bond, it is removed from the "active list." A new spin is added 
to the list as soon as it acquires a broken bond due to a transition of one of its neighbors. The memory 
required for the list is proportional to the length of the interface, L x /\cos(f>\. No lattice boundaries are 
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needed in the y-direction, and consequently arbitrarily long simulation runs can be performed. Except for 
these modifications, the algorithm is a straightforward implementation of the discrete-time n-fold way 111]]. 

The mean tilt angle <p w & s fixed by helical boundary conditions in the cc-direction. The production runs 
were performed with L x = 1000 and fixed <fi between and 7r/4. Since kinetic Ising interfaces belong to 
the KPZ universality class || ||, the macroscopic interface width, ([y(x, t) — (y(x 1 t))} 2 ) 1 ^ 2 , is expected to 
saturate at a value proportional to L]J 2 after a time t oc L% . In order to ensure stationarity we ran 
the simulation for 10 000 MCSS before taking any measurements. Exploratory runs with both larger and 
smaller L x and "warm-up" times showed that the values used in the production runs were sufficient to 
ensure a macroscopically stationary interface. Class populations and interface velocities were averaged over 
200 000 MCSS. The macroscopic stationarity of the interface should abundantly ensure the stationarity of 
these quantities, which are properties of the local interface structure. 

5.2 Interface velocities 

The overall quality of the four different approximations developed above (linear-response (Sec. ^|) and 
nonlinear-response (Sec. ||), each either excluding or including the SOS- violating transitions from classes 
10s and 20s) were explored in a wide range of H and T. In Fig. |^ the four approximations are compared 
with MC data for interfaces parallel to the x-axis at T = 0.2T C [Fig. |(a)] and T = 0.6T C [Fig. |(b)], 
where T c is the exact critical temperature. The crossover field Hmfsp(T) is marked by a vertical dashed 
line in both (a) and (b). At the lowest temperature, the interface has very few steps higher than one for 
H = 0. As a result, it makes practically no difference whether SOS-violating transitions are allowed or 
not, and the curves representing the linear-response approximations including and excluding SOS-violating 
transitions coincide in Fig. ||(a). The nonzero field increases the average step height, and as a result the 
SOS- violating transitions contribute significantly for stronger fields in the nonlinear-response approximation. 
Overall, the nonlinear- response approximation including SOS-violating transitions ("nonlinear inclusive" or 
NLI) is everywhere better than the others and is particularly good for H < Hmfsp(T). It is the only one 
which will be used in the rest of this paper. 

The temperature dependence between T = and T c of the velocities of interfaces parallel to the x-axis 
are shown in Fig. || for several values of H/ J between 0.2 and 2.0. For H/J < 1.5, the discrepancy is nowhere 
greater than a few percent - mostly much smaller. Even for H/J — 1.9 and 2.0, the discrepancy remains 
below about 15% everywhere. 

The anisotropy of the interface velocities is shown in Fig. ||for several values of H/J between 0.1 and 2.0. 
The agreement between the NLI approximation and the simulations is very good, except for the strongest 
fields. At T = 0.2T C [Fig. ||(a)], both the analytic and simulation results for H/J < 1.5 increase with 
increasing <f>. As predicted by Eqs. ([uj) and (jl^), under weak fields they cross over from the form of the 
PNG model for tan0 -C X, to almost exact agreement with the single-step model for tan0 3> X [Fig. ^(b)]. 
Growth shapes generated from these velocities would be almost square, with their sides parallel to the x and 
y axes 0. For strong fields, however, the analytic approximation predicts velocities that are slightly larger 
in the symmetry directions than for inclined interfaces, as is the case for the Eden model |36|, |37], This 
reverse anisotropy is not seen in the MC data, which become almost perfectly isotropic and would lead to 
circular growth shapes. With the Glauber dynamic it appears that stronger fields and lower temperatures are 
needed to observe reverse anisotropy in the MC data ||. However, the discrepancies between the simulations 
and the analytic approximation arc modest in the regime explored here, even near <j> = 7r/4 for strong fields. 
For small <fi the agreement remains good for all fields H/J < 1.9. At T = 0.6T C [Fig. g(c)], the simulated 
velocities are practically isotropic for all fields. The analytic NLI approximation works well, except for 
gradually increasing Eden-like reverse anisotropy for the stronger fields. 

5.3 Spin-class populations 

While our main emphasis is on estimates of the interface velocity, the most detailed information about 
the strengths and weaknesses of our mean-field approximation for the interface structure is found in the 



9 



individual spin-class populations. Examples of the ii-dependences for cf> = at T — 0.6T C are shown in 
Fig. 0. Individual populations are shown in Figs. 0(a) and (b), while spin classes with the same number 
of broken bonds are combined in Fig. 0(c). The crossover field -ffMFSp(r) is shown in all three panels as 
a vertical dashed line. The agreement between the theoretical estimates and the MC data is good out to 
about this field. 

The discrepancies that develop with increasing H stem from two sources. One is the fact that the 
simulated interface, unlike the SOS description used in the theoretical analysis, is not restricted to be free 
of overhangs and bubbles. The other is the development of correlations between neighboring step heights, 
which are not included in our single-step mean-field theory for the local structure of the driven SOS interface. 

The presence of overhangs and bubbles is clearly reflected in the increase of (n(01s)), relative to the 
monotonically decreasing analytic approximation [Fig. ujfa)], and in the nonzero populations of the classes 
that are not populated in the SOS approximation [Fig. 0(b)]. 

The increasing correlations between nearest-neighbor steps is expressed by the gradual disappearance of 
the symmetry between the populations of spins with s = +1 (stable) and s = — 1 (unstable) as the field 
is increased. The effect of these correlations is typically to broaden protrusions on the leading side of the 
interface ( "hilltops" ) and sharpen protrusions on the trailing side ( "valley bottoms" ) |}2j , or vice versa |3j| . 
In terms of spin-class populations, the former corresponds to (n(21+)) < (n(21— )) and (n(ll+)) > (n(ll— )). 
This is precisely what is seen in the simulations [Fig. 0(a)]. More generally, we find that the populations of 
s = +1 arc enhanced for the classes with one and two broken bonds, while the populations of s = —1 are 
enhanced for the classes with three and four broken bonds [Fig. 0(c)]. The latter include both sharp valley 
bottoms and bubbles of the unstable phase, which form a wake behind the moving interface. This shows 
that the breakdown of the SOS description and the evolution of lateral correlations are not independent at 
strong fields. The overall picture is essentially the same at T = 0.2T C : both the populations in the non-SOS 
classes, and the asymmetry become significant only near Hmfsp(T), which at that temperature is close to 
H/J = 2. 

6 Discussion 

In this paper we introduced and explored analytic approximations for the propagation velocity and spin-class 
populations of a field-driven interface in a two-dimensional kinetic Ising model. The model is applicable to 
a number of systems undergoing order-disorder phase transformations, including magnetic and ferroelectric 
systems and crystal growth under conditions which are not limited by diffusion. 

The approximations are based on a linear-response approach, in which the equilibrium fluctuations of 
the interface (as embodied in average spin-class populations) were estimated from the Burton-Cabrera- Frank 
SOS model. The theory was extended by developing a mean-field approximation for the field-dependent spin- 
class populations in a moving flat interface, yielding a nonlinear-response approximation for the interface 
velocity. This extension considerably improves the agreement with MC simulations. 

Our simulation results are consistent with those of Ref. Q. However, since that study used the Metropolis 
dynamic [see their Eq. (3)], the velocities cannot be compared directly. For asymptotically low T, the velocity 
is completely determined by the kink density, (n(lls)). This is the regime in which the single-step and PNG 
models are expected to hold. Indeed, they emerge in the proper limits from our analytic approximation, as 
shown in Eqs. ([l4|) and (|l5|). For larger H <2J and T < T c , other spin classes also contribute significantly 
to the interface velocity. In this parameter range we obtain good agreement between theory and simulation 
almost everywhere. This agreement includes the disappearance of the velocity anisotropy with increasing T 
and H. 

Both the SOS approximation for the interface structure, and the assumption of uncorrelated step 
heights employed in the nonlinear-response approximation break down for stronger fields. While the 
NLI approximation for the interface velocities remains very satisfactory overall, the breakdown of these 
assumptions can be detected more clearly in the spin-class populations. The only detailed theories of 
the local structure of driven interfaces that we are aware of, are for much simpler SOS models than the 
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unrestricted one studied here, such as the body-centered SOS (BCSOS) and restricted SOS (RSOS) models 
[ p2[ . It would therefore be interesting to compare our mean- field approximation for the interface structure 
with MC simulations of the driven SOS interface, so that the complications arising from overhangs and 
bubbles can be avoided Q . 

In summary, the approximations presented here give results for the spin-class populations and propagation 
velocities of Ising interfaces driven by nonzero fields, which are exact in the asymptotic limits of low T and 
H, and which agree well with MC simulations almost everywhere, even for H near 2J and T near T c . It is 
not clear whether the wide regime of applicability is an accident limited to two dimensions, and it would 
therefore be interesting to extend the approach to three dimensions as well. 
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Table 1: The spin classes in the anisotropic square-lattice Ising model. The first column gives the class 
labels, jks. The second column gives the total field and interaction energy for a spin in each class, E(jks), 
relative to the energy of the state with all spins parallel and H — 0, -Bo = —2(J X + Jy)- The third column 
gives the change in the total system energy that would result from reversing a spin in a particular class 
from s to — s, AE(jks). In both E(jks) — Eq and AE(jks), the upper sign corresponds to s — —1, and the 
lower sign corresponds to s = +1. The first three classes (marked *) have nonzero populations in the SOS 
approximation, and flipping a spin in any of them preserves the SOS interface configuration. The next two 
classes (marked t) also have nonzero populations in the SOS approximation, but flipping a spin in any of 
them may produce an overhang or a bubble. The two classes marked \ have zero populations in the SOS 
approximation, but flipping a spin in any of them may lead to a configuration compatible with the SOS 
constraint. Class 22s represents a single spin which is antiparallel to all its neighbors; flipping such a spin 
yields a bulk spin in class 00— s. Although only the classes marked * and f have nonzero populations in 
the SOS approximation, the MC transition probabilities of all classes except 00s are given by Eq. (Q). The 
bulk spins, 00s, have zero transition probabilities in the dynamic used here. 



Class, jks 


E(jks) - Eo 


AE(jks) 


01s * 


±H + 2Jy 


T2H + AJ X 


lis * 


±H + 2(J X + Jy) 


T2H 


21s * 


±H + 2(2J X + J y ) 


T2H - AJ X 


10s f 


±H + 2J X 


T2H + AJy 


20s f 


±H + 4 J x 


T2H - A(J X - Jy) 


12s % 


±H + 2(J x + 2J y ) 


T2H - AJy 


02s X 


±H + 4 Jy 


T2H + A(J X - Jy) 


22s 


±H + 4(J X + Jy) 


T2H - 4(J X + Jy) 


00s 


±H 


T2H + 4{J X + J y ) 
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Table 2: The mean populations for the spin classes included in the SOS approximation, together with 
their contributions to the interface velocity under the Glauber dynamic. The first column gives the class 
labels, jks. The second column gives the mean spin-class populations for general tilt angle 4>, with cosh 7(0) 
from Eq. (^|). The third column gives the spin-class populations for (f> = 0. Using X = e~ 2,3Jx in these 
expressions, one obtains the linear-response result in which the spin-class populations are evaluated for 
H = 0. Using X = X(T, H) from Eq. (Efl) with the transition probabilities of the particular dynamic used 
[here: Glauber, Eq. (Ell)], one gets the nonlinear-response approximation, which includes an estimate of 
the field-dependent modifications of the spin-class populations in the moving interface. The fourth column 
contains the contributions to the mean interface velocity in the y-direction from spins in classes jk— and 
jk+, using the Glauber dynamic. These are the only quantities in this table that depend explicitly on the 
specific dynamic, and they could easily be replaced with results for, e.g., Metropolis. 



Class, jks 


(n(jks)), general <fi 


(n(jfcs)), </> = 




01s * 


l-2Xcosh7(0)+A"^ 




tanh(/3i?) 




(i+x) 2 


, , rainh(23Jx)l 2 
1+ co S h((JH) 


lis * 


2X[(1+X 2 )cosli7(0)-2X] 
(1-X 2 ) 2 


2X 
(1+X) 2 


tanh (j3H) 


21s * 


X 2 [1-2X cosh~/(<f>)+X 2 ] 


X 2 


tanh{f3H) 


(l-X 2 ) 2 


(1+X) 2 


1 + 


"sinh(2f).7 :c )l 
c(»sli{ iiH) 


2 


10s f 


2X 2 \ 2cosh 2 7(0)-l-2Xcosh7(4>)+X 2 


2X 2 (1+2X) 


tanh(/3if) 


1-X' 2 [ l-2Xcosh7(0)+X 2 

X 2 [1-2X cosh7(0)+X 2 ] \ 
(1-X 2 ) 2 J 


(1-X 2 )(1+X) 2 


1 + 


sinh(20J y ) 
co 3 h(/3J7) 


2 


20s f 


X 4 [1-2X cosh7(0)+X 2 ] 


X i 


tanh(/3ff) 


(1-X 2 ) a 


(1_X 2 )(1+X) 2 


1 + 


Binh(2^(J x -Jy)) 1 2 
cosh(f3H) 
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Figure 1: A short segment of a zero-field equilibrium SOS interface between a positively magnetized phase 
for y < and a negative phase for y > 0. The independent step heights, 5(x), are drawn from the pdf given 
by Eqs. Q and (^) with T = 0.6T C and <j> — 0. Interface sites representative of the different spin classes 
compatible with the SOS approximation are marked with the notation jks explained in Sec. ^. Sites in the 
uniform bulk phases are 00— and 00+. 
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Figure 2: Comparison between MC simulations and the four different analytic approximations introduced 
here. Shown is the normal interface velocity (v±) vs H for <j> = at T = 0.2T C (a) and T = 0.6T C (b). MC 
data (squares), linear-response excluding SOS-violating transitions (dotted), linear-response including SOS- 
violating transitions (short-dashed), nonlinear-response excluding SOS-violating transitions (long-dashed), 
and nonlinear-response including SOS-violating transitions (NLI) (solid). In the low-temperature case 
shown in (a) the two linear-response approximations are almost indistinguishable, and the dotted and 
short-dashed curves coincide as the lowest of the three curves shown. Statistical errors in the MC data 
are everywhere much smaller than the symbol size, both in this and all subsequent figures. The vertical 
long-dashed lines mark the crossover field Hmfsp(T). In both (a) and (b), the NLI approximation gives 
the best overall agreement with the data. It is the only one which will be used in subsequent figures. 
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Figure 3: Figure for calculating the transition rates W[5(0) -> 5(0) + 1] and W[5(0) + 1 -> 5(0)], Eq. @. 
Interface configurations are indicated by bold line segments. Like in Fig. [l], the negatively magnetized phase 
is above the interface, and the positively magnetized phase below it. At the center is shown a step 5(0) > +1 
(here shown as 5(0) = +1 for concreteness). A transition to (5(0) + 1) can be effected by flipping either one 
of the spins in the dashed boxes, each with probability 1/2. Flipping the initially negative spin (— ) decreases 
5(+l) by one. The resulting configurations are shown to the right. The energy change depends on the initial 
value of 5(+l), and the two possible energy changes and the corresponding conditions are shown next to 
the right-pointing arrows. The thin horizontal lines represent the interface configuration corresponding to 
the equality in the conditions. The energy changes and conditions for the reverse transitions are shown 
next to the left-pointing arrows. Flipping the initially positive spin (+) analogously decreases 5( — 1) by 
one. The energy changes and conditions for the forward and reverse transitions resulting from flipping this 
spin, are indicated in the left-hand half of the figure. 
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Figure 5: Temperature dependence of the normal interface velocity, (v±), for (j> — at H/J = 2.0, 1.9, 1.5, 
1.0, 0.5, and 0.2 from top to bottom. MC simulations (data points) and the NLI analytic approximation 
(solid curves). 
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Figure 6: The normal interface velocity (v±), shown vs tan <p for different fields. Simulation data are shown 
as data points and the NLI approximation as solid curves. The fields included are (from top to bottom) 
H/J — 2.0 (empty squares, only in a and c), 1.9 (empty triangles, only in a), 1.5 (filled diamonds), 1.0 
(filled stars), 0.5 (filled triangles), and 0.1 (filled squares), (a): At T = 0.2T C . (b): The four weakest 
fields at T = 0.2T C , divided by tanh(/3i?). This shows the crossover between PNG and single-step growth, 
Eqs. (0) and @. The dotted curve is the zero-temperature single-step result, (c): At T = 0.6T C . 
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Figure 7: Comparison between MC simulations (filled symbols for s = +1 and empty symbols for s = — 1) 
and the NLI analytic approximation (solid curves) for the spin-class populations. Results for T = 0.6T C 
and <f> = are shown vs H. (a): Individual populations of the spin classes compatible with the SOS 
picture. From top to bottom in the left part of the figure are 01s (stars), lis (diamonds), 10s (stars), 21s 
(triangles), and 20s (diamonds), (b): Individual populations of the spin classes that would be unpopulated 
for a perfect SOS surface: 02s (diamonds), 12s (triangles), and 22s (squares). These nonzero populations 
clearly indicate the gradual breakdown of the SOS approximation with increasing H. (c): Class populations 
combined according to number of broken bonds, j + k. From top to bottom: one (stars), two (diamonds), 
three (triangles), and four (squares) broken bonds. 
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